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Abstract 

A protein model with the pairwise interaction energies varying as local en- 
vironment changes, i.e., including some kinds of collective effect between the 
contacts, is proposed. Lattice Monte Carlo simulations on the thermody- 
namical characteristics and free energy profile show a well-defined two-state 
behavior and cooperativity of folding for such a model. As a comparison, 
related simulations for the usual Go model, where the interaction energies 
are independent of the local conformations, are also made. Our results in- 
dicate that the evolution of interactions during the folding process plays an 

important role in the two-state cooperativity in protein folding. 
PACS number: 87.15 Aa, 87.15 Cc, 87.15 He 
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Due to the developments of experimental methods and theoretical models, many achieve- 
ments of protein folding have been made recently . A protein can fold itself to its uniquely 
well-defined native structure in a biologically short time, regardless of the huge number of 
possible conformations, showing a highly cooperatively kinetic behavior. It is now clear that 
the cooperativity of folding may result from the backbone hydrogen bonding, sidechain pack- 
ing and hydrophobic interactions, among them the hydrophobic interactions are believed to 
be the dominant driving force for folding For many small single-domain proteins or 
lattice proteinlike models, there is a two-state behavior between the unfolded states and the 
folded native one PJ^. Recently, Chan and Kaya indicated that according to the calori- 
metric criterion, which is widely used in experiments as a condition for two-state folding, 
popular lattice models, e.g., two-letter HP and 20- letter MJ models, are far from two-state 
models. This may be due to some flawed assumptions in the potential functions used in these 
models. Lattice models usually use statistical potential functions extracted from the pair- 
ing frequencies of 20 kinds of amino acids in databases of protein structures 0. Although 
these knowledge-based potentials may be a good approximation to the relative strength of 
interactions between the residues in the native state, they provide no information about 
how the interactions evolve during the folding. For computational convenience, a common 
assumption in lattice models is that the interactions are additive, and they are the same 
during the folding as in the native state. This means that the interaction energies are 
conformation- independent. Clearly this is not relevant to the experimental situation 0. In 
fact, as Dill pointed out ||^, the thermodynamic additivity principle which is widely used in 
chemistry may be unsuitable in biochemistry. Some recent experiments also indicated that 
the transition state is an expanded version of the native state, where the majority of inter- 
actions are partially formed 0, and their strengths are different from those in the native 
state (with $ < 1). That is, these interactions depend on the conformations [0, especially 
the local structures around the contacts as emphasized recently in Ref. |T0[. Previously, the 
non-additivity was built in a lattice model for packing effects [[T^]; and the hydrophobic force 
depending on the local density of peptide atoms was also taken into account in an off-lattice 



model |T^. Studies on these models show that the introduction of the non-additivity is 



significiant, but the two-state cooperativity of these models is not checked, and the effects 
of the non-additivity on thermodynamics and kinetics of folding need to be further studied. 

In this paper, we develop a refined Go model where the pairwise interaction energies 
vary as the local environment changes, i.e., some kinds of collective effects between contacts 
are introduced. Our purpose here is to study the two-state cooperativity of protein folding 
and its physical origin with such a model. Our results give a general picture about how the 
conformation-dependent interactions affect the folding kinetics, which is consistent with the 
phenomenological explanation based on experimental results. 

We model a polypeptide chain as a self-avoiding chain on a cubic lattice. A contact is 
formed if two residues are space adjacent but not sequence adjacent. If two residues form a 
contact as the same as in the native state, we call this contact a native contact, otherwise 
a non-native contact. Following the Go model |T^, only native contacts are considered to 



contribute to the total energy. Different from the Go model, we assume that the interaction 
energies between residues are conformation-dependent, and vary with changes of the local 
environment. To achieve this, we introduce a parameter S to describe the degree for a residue 
being ordered relative to the native state. For the i-th residue in a certain conformation, its 
degree of order Si is defined as 

S^ = Z,IZT\ (1) 

where Zi is the number of native contacts in this conformation, 2;""* is the number of contacts 
formed in the native state. Obviously, Si varies between (the i-th residue being fully 
disordered) and 1 (being fully ordered). Thus, the interaction energy between residues i and 
j, Bij=-{Si+Sj)e /2 is defined, where Si and 5*^ are the degrees of order for residues i and j, 
respectively, e is the unit of energy and is set to be 1 in this work. The total energy of the 
conformation then is E = Yl,i<j ^ijBij, where A is unity when residues i and j form a native 
contact, and zero otherwise. Here, a contact formed between residues i and j may have 
different energies in different conformations, i.e., Bij may change from one conformation to 



another (for the Go model, one always has Bij = —e). In general, a contact formed between 
residues i and j will stabilize, to some extent, other contacts that residue i or j formed 
with other residues. On the contrary, its breakage may destabilize those contacts as well. 
Therefore, the introduction of the degree of order for a residue into the potential function 
reflects the cooperativity between the residues. Although the correlation distance is small, 
only one lattice unit, the many-body effects are obviously included in our model. Figure 1 
shows such a collective effect. The interaction energies of contact A-B (or B-C) are different 
when the other contact is present or not present. Clearly, the energy of state I3 is lower than 
the sum of that of states Ji and I2, indicating the interaction non-additivity. Each contact 
is stabilized by the other contact due to the collective effect. Note that in this paper our 
model is called Go+ model to distinguish from the Go model. 

Now let us present the Monte Carlo simulations on the thermodynamic and kinetic 
features for both models. The mean first passage time (MFPT), as a common measure of 
folding rate, is calculated by an average of the first passage time (FPT) over 1000 runs. 
Each run begins with a random conformation, and ends when the native state is reached for 
the first time. The FPT is the Monte Carlo steps (MCS) consumed in a run. 

Generally, as the temperature T decreases, the population of the native state, P/v, in- 
creases from zero to about unity. The degree of sharpness of changes in P/v, similar to the 
"rapidity" in Ref. fl^, is a measure of the cooperativity of the folding reaction. Figure 2 



shows the population P/v and the specific heat Cy versus temperature T for a 36-mer chain 
for both models. P/v is defined as P/v = e~^^ /'^ / ^E^{E)e^^^'^ ^ where Vl{E) is the density 
of states for energy E, En is the energy of the native state. Q{E) is calculated with the 
Monte Carlo histogram method [|15[. From Fig. 2 we can see that the folding transition for 



our Go+ model is much sharper than that of the Go model, i.e., a sharper change in P/v. 
There is also a single peak in curve, but it is narrower than that of the Go model. For 
our Go+ model, the maximum of occurs at a temperature nearly the midpoint tempera- 
ture of transition with P/v=l/2, i.e., the difference between these two temperatures is quite 



small. This is consistent with recent studies on naturally occuring proteins |1J,|16[, implying 



a good cooperativity of folding in this model. Differently, such a temperature difference is 
large for the Go model (see Fig.2), indicating that the folding of the Go model is much less 
cooperative than that of the Go+ model. Since the sharpness is only a qualitative descrip- 
tion for the transition, we further calculate the equilibrium energy distribution at the folding 
transition temperature, Tf. Figure 3 shows such distributions for both models. Clearly our 
Go+ model shows a good bimodal behavior, and the denatured-state energy is distributed in 
a narrow region [see Fig. 3(a)]. This means clearly a two-state folding and there is basically 
no intermediate states at equilibrium. Differently, for the Go model as shown in Fig. 3(b), 
there are many intermediate states and the bimodal behavior is not so significant as that in 
Fig. 3 (a). Thus for the Go model the folding is not of a two-state. This is in agreement with 
Chan and Kaya's argument 0. 

In experiments, a well-established criterion for two-state folding is that the van't Hoff 
enthalpy AH^h around the transition midpoint is equal, or very close, to the calorimetric 
enthalpy AHcai of the entire transition. In this work, we calculate the ratio AH^h / AH^ai as 
suggested in Ref. |^ (here, the definition of AH^h / AHcai is equal to (A;2)^ in Ref. |^), and 
list the results in Table I. From Table I, we can clearly see the difference between the Go 
model and our Go-|- model. The Go model, which is considered as a model with minimal 
energetic frustrations, does not meet the calorimetric two-state criterion and gives out the 
value of AHyn/AHcai far from 1. Nevertheless, our model satisfies the criterion quite well 
(for real proteins, the value of AH^h / AH^ai is 0.96±0.03 fl^])- This, again, implies the 
two-state folding and the good cooperativity of our Go-|- model. 

Physically, the high cooperativity of our model may result from the narrow distribution of 
the denatured states and the high population of the native state at the folding temperature 
(see also Figs. 2 and 3). In our model, the energy spectrum relating to various conformations 
is redistributed, comparing with that of the Go model, due to the collective effect between 
interactions. As a result, the energies of non-native conformations are moved to higher 
energy levels and a larger energy gap is left between the non-native conformations and the 
native one (for the two models, the energies of the native state are the same). The large 

5 



energy gap makes the native state paticularly stable, which is beheved to be a necessary 
condition for cooperative folding [|T^. This may be the physical origin of the two-state 
cooperativity. It can be further explained from the viewpoint of the free energy profile. For 
our Go+ model, as shown in Fig. 4, the free energy profiles have broad activation barriers. 
The broad activation barriers can account for the large movement of transition state caused 
by mutation or temperature changes, and are considered as a common feature of the two- 
state folding . Our numerical results are surprisingly consistent with a phenomenological 
speculation for the existence of such a free energy profile in Ref . . It should be noted that 



the broad activation barriers are consistent with the narrow distribution of the denatured 
states. 

Now let us make a comparison of the foldability based on the plots of the MFPT versus 
Pn for both models. Note that we use P/v instead of the commonly used temperature T 
in the horizontal axis in Fig. 5. This is because that an identical condition should be taken 
for the comparison. In lattice simulations, the temperature has an arbitrary unit and also 
has no direct relationship with the real temperature. The comparison between two different 
models at the same temperature may make no sense. Nevertheless, at an identical condition 
with the same P/v, the differences in the foldability can be well-defined. This is similar to 
other conditions used previously From Fig. 5, we can see that the MFPT for our Go+ 
model shows a slow decrease as Pjv increases, it reaches a minimum at Pjv ~ 0.93, and then 
it increases. For the Go model, there is also a minimum but at P/v ~ 0.71. It is clearly 
that when the native state is stable (say, P/v > 0.9), our Go+ model folds significantly fast, 
i.e., the MFPT is smaller with one or two orders of magnitude than that of the Go model. 
Physically, this can be explained as follows. From Eq.(l) we can easily see that the energy 
gain of forming a contact is usually smaller for our Go+ model than that for the Go model. 
At high temperatures, entropic contribution is dominant to the free energy barrier, and 
the loss of entropy is always undercompensated by the energy gain, thus the Go-|- model 
folds slower for its smaller energy gain. Whereas at low temperatures, folding is nearly a 
downhill process, and the loss of entropy is always overcompensated by the energy gain. 
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Therefore, for the Go+ model, it is easier to escape from kinetic traps, and the folding is 
faster. Finally, we note that for the two models the pathways of reaching the transition state 
from the denatured state are different. Due to the high cooperativity in our Go+ model, a 
good core, the assembly of non-polar residues, is formed much earlier at low temperatures 
than that in the Go model. Detailed kinetic results will be reported elsewhere. We also note 
that similar results are obtained for different chain sizes. 

In conclusion, our Go+ model, with many-body interactions depending on the local 
structures included, exhibits a good two-state folding behavior. Our results suggest that 
the evolution of interactions during the folding plays an important role in the two-state 
cooperativity in protein folding. We give a possible way how the interactions evolve in 
the folding, which may capture some essential features of the two-state folding. We expect 
further study could provide new insights into the mechanism of protein folding. 

We thank H.S. Chan, A. Maritan and D. Thirumalai for useful suggestions. This work 
was supported by the Foundation of NNSF (No. 19625409, and No. 10074030). 
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Table I: The ratios of AH-^h / ^Hcai for the Go model and our Go+ model, respectively. 
Ten sequences are calculated for each chain size. 
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FIG. 1: Schematic illustration of collective effect between two interactions. From a state 
Jo with three unstructured residues, the chain can be settled in a state Ji (or I2) with a 
contact A-B (or B-C) and an equilibrium constant Ki (or A'2). A state Is with two contacts 
A-B and B-C can be reached from state Ji or I2, but with different equilibrium constants 
i^27 or In state I3, each interaction is stronger by a factor 7 due to the existence of 

the other contact. 

FIG. 2: Population Pjv and specific heat varying with the temperature T for a 36-mer 
chain. 

FIG. 3: The energy distribution for the same 36-mer used in Fig. 2, using (a) Go+ potential 

and (b) Go potential at respective folding transition temperature, Tf. 

FIG. 4: The free energy profile F{E) = E — TS{E) of our Go+ model at different 

temperatures, where entropy S{E) is calculated by using entropy sampling Monte Carlo 

method |21|]. Here U, N and TS denote the unfolded state, native state and transition state, 

respectively. Note that the free energy profile at high temperature is overall shifted so that 

the unfolded states are overlapped. 

FIG. 5: MFPT versus P/v for a 36-mer chain. 
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. . AH,H/AHeal 

Cham size 



Go model Go+ model 

27-mer 0.62+0.01 0.91+0.01 

36-mer 0.60+0.01 0.91+0.01 

48-mer 0.74±0.01 0.95±0.01 



